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A new approach to account for finite thermal conductivity and turbulence effects within 
atomizing droplets of an evaporating spray is presented in this paper. The model is an extension of 
the T-blob and T-TAB atomization/spray model of Trinh and Chen [9]. This finite conductivity 
model is based on the two-temperature film theory in which the turbulence characteristics of the 
droplet are used to estimate the effective thermal diffusivity for the liquid-side film thickness. Both 
one-way and two-way coupled calculations were performed to investigate the performance of this 
model against the published experimental data. 

Introduction 

Spray vaporization and combustion studies are of primary importance in the prediction and 
improvement of systems utilizing spray injection. In liquid fuelled combustion systems such as 
industrial boilers, gas turbines, direct ignition diesel engines, rocket and air-breathing engine 
applications, the combustion performance is highly dependent on effective liquid fuel atomization 
and its subsequent evaporation processes. The recent review by Sazhin [1] had summarized the 
development of various droplet heat-up and evaporation models with different levels of complexity. 
As also pointed out in [1], to numerically simulate evaporating spray using CFD (computational 
fluid dynamics) methodologies, detailed droplet heating and evaporation models based on single- 
droplet analysis have to be simplified to be CPU efficient. Although the film-theory based 
evaporation models had been extensively used, most of the developments were done on the gas-side 
heat and mass transfer [2, 3]. Within the liquid droplet, the two-temperature formulations, in which 
the finite-conductivity (F-C) effect was modelled by the temperature difference between the droplet 
surface temperature and a droplet “core” temperature have been proposed recently. In the model of 
Renksizbulut et al. [4], the difference between the surface and core temperature was related to the 
heat flux at the droplet surface by a constant Nusselt number. Zeng and Lee [5] developed a zero- 
dimensional model, in which the difference between surface and core temperature was traced by an 
ordinary differential equation (ODE) to account for the non-uniform distribution of temperature 
inside the droplet. The ODE two-temperature model was also used by Miller et al. [6] to account for 



non-equilibrium Langmuir-Knudsen mass transfer effect. Ra and Reitz [7] used a thermal boundary 
layer within the droplet to account for the finite-conductivity effect. The thermal boundary layer 
thickness was calculated using the thermal diffusivity model of Abramzon and Sirignano [8] to 
account for the droplet internal circulation. 

In this paper, a new finite-conductivity model is developed upon the two-temperature 
formulation. The finite-conductivity effect is phenomenologically modelled through a thermal 
boundary layer within the droplet. The thermal diffusivity is calculated based on the turbulent 
characteristics within the droplet. The current study is an extension of a recently developed 
atomization/spray model [9, 10], the T-blob/T-TAB model, to include spray evaporation effects. 

J 

Due to the unique feature of T-blob/T-TAB, in which the turbulence characteristics is accounted for 
within the droplet phase, extension of this model to include finite conductivity effect in the 
evaporating droplet can be made naturally. The model development will be described in this paper. 
Validations for one-way one-dimensional and two-way coupled multi-dimensional evaporating 
sprays will be presented. 

Model Development 

The current vaporization model is developed for computational analysis based on the 
Eulerian-Lagrangian numerical approach. In this formulation, the spray/droplets dynamics is 
described in a Lagrangian coordinate such that numerical droplets are tracked within the Eulerian 
gas dynamics. Liquid phase is tracked from the injector plane, and the primary atomization as well 
as the subsequent secondary break-up is modelled using the T-blob/T-TAB hybrid model of Trinh 
and Chen [9]. Both primary and secondary droplet break-up processes are modelled and the 
transition from primary to secondary break-up is modelled based on an energy balance. In addition 
to the droplets position and velocity, liquid turbulence is accounted for through the injector 
characteristics by the two-equation k-e turbulence model formulation using the T-blob/T-TAB 
model. The inherent turbulence in the injected fuel spray affects the heat and mass transfer rates of 
the vaporization process. The effects of these changes in the rates have to be accounted for in the 



numerical models for spray evaporation. Detailed model description and validations can be found in 
[9, 10], and it is sufficient to say that within each numerical droplet, turbulence characteristics such 
as fluctuating velocity level, length and time scales are supplied by the model. 

To utilize the T-blob/T-TAB model liquid j et atomization, turbulence characteristics need to 
be specified as the inlet boundary conditions. Based on integral analysis of straight injector [9], 
liquid turbulent kinetic energy and its dissipation rate at the injector nozzle exit are estimated from: 



where L is the length of the injector nozzle, D nozz ] e is the nozzle diameter, and the jet injection 
velocity U. A set of ODE’s were derived to track the evolution of ki and £/ within the droplet 
according to the T-blob/T-TAB model. The values obtained from the evolution of ki and €/ are used 
in the heat transfer calculations of the evaporation droplet. 

To relax the infinite-conductivity assumption, and thus perfect mixing within the liquid 
droplet, a “two-temperature” model is formulated. In the two-temperature model, the core (or bulk) 
temperature (Td), is assumed to be well-mixed by convection/turbulence transport. Consistent with 
the “film theory”, heat resistance exists at the near surface region, and the droplet surface 
temperature (T s ) differs from the droplet core temperature. The heat transfer coefficient across this 
thin film (or boundary layer) is then formulated through the turbulence characteristics supplied from 
the T-blob/T-TAB model, to account for the finite conductivity effect. In the Lagrangian coordinate, 
the heat-up of the droplet core is formulated as: 

m d c p ^= h ,{T s -T d ) A * ( 3 ) 

where, hj is the liquid-side heat transfer coefficient and Ad is the droplet surface area. The heat 
transfer coefficient is determined from the thermal conductivity and a thermal boundary layer 


formulation [11] as: 



h, =§*-, ( 4 ) 

°e 

where \i is the liquid thermal conductivity. Using an unsteady equivalent boundary layer thickness, 
the film thickness, o e , is given by a cff t . The time scale t is estimated based on the heat transfer- 

limited integration time step, and should be independent of the global time step used in numerical 
calculations (see the results in figure 1). In this study, we adopted the evaporation sub-cycle time 
step based on the formulation described in Amsden et al. [12], This time step is computed by the 
relation given by: 


P g M g (4xr d )N p 


( 5 ) 


where dV is the volume of the computational cell, p is density, p is the viscosity of gas mixture 


within the cell, and N is the number of droplets for each computational particle p. The effective 

thermal diffusivity (ofeff), based on the turbulence characteristics within the droplet, is estimated 
from: Qfeff = G!i am + QWb, in which the turbulent thermal diffusivity is calculated from the two- 


equation turbulence model diffusivity formulation: 


a ,urb = 


C k 2 

Z£ K 1 

Pr, s. 


(6) 


here, Pr, is the turbulent Prandtl number and is set to be 0.9. The liquid droplet turbulence 
quantities k, and s, are obtained from the T-blob/T-TAB atomization/spray model. In cases where 

gas phase diffusivities are much larger than liquid diffusivities, the droplet core heating will be rate 
controlling and the gas-side heat/mass transfer will respond in a quasi-steady manner. The surface 
temperature of the droplet is determined from a heat and mass transfer balance at the interface 
between the droplet and the surrounding gas assuming no heat accumulation at the droplet surface 
such that: 


K m d = 9,-9, 


( 7 ) 



where L v is the latent heat of the fuel at the surface temperature, and qg is the heat transfer rate from 
the environmental gas to the surface. In this paper, the classical Spalding evaporation model is used 
to model the gas-phase transport, thus the gas heat transfer rate was calculated as: 


q g = 2nr d A g Nu g (T g -T s ) 


ln(l + 5 r ) 

B t 


( 8 ) 


qi in Eq. (7) is the heat transfer rate from the droplet interior to droplet surface, and is calculated 
from: 


-MLiIA 


( 9 ) 


t he evaporation rate at the droplet surface is given as: 

- = 27rr, (p g D)Sh g !n(l + BJ (10) 

In the above equation, B m is the Spalding mass transfer number, Sh g is the Sherwood number, and 
D is the binary diffusivity. The Sherwood and the Nusselt numbers were calculated using the 
classical correlations give by: 


Nn =2 + 0.6 Re 05 Pr 0JW 

S P 


( 11 ) 


Sh g = 2 +0.6Re p 05 Sc om 

( 12 ) 

The solution algorithm used in this study starts with an estimated surface temperature (T s ) at a new 
time step. The Clausius-Clayperon equation and the Raoult’s law are then used to calculate the fuel 
vapour molar fraction followed by the calculation of the evaporation rate. An estimation of Ta is 
also required to simultaneously satisfy Eqs. (3) and (7). A more detailed iterative procedure can be 
found in Balasubramanyam [13]. 

Results and Discussion 

One-Way Evaporating Atomizing Spray 

In one-way simulation, the gas flow field was prescribed and the droplets were tracked in time 
domain. This case is an extension of the one-way T-blob/T-TAB testing case described in Trinh and 



Chen [9]. Iso-octane fuel similar to the fuel used in reference [7] was issued through a long injector 
tube at 300°K. The length of the injector nozzle (L) is 0,8 mm and the nozzle diameter (D n0Z zie) is 
0.3 mm. A jet injection velocity (U) of 102 m/s was used for the test case. The environment is 
quiescent nitrogen at a temperature of 600°K. The gas properties were calculated based oh the 
reference state determined from the ‘l/3 rd ’ rule [11], In this calculation, a ‘blob’ of numerical 
droplet was injected at the orifice plane with orifice diameter. The droplet then went through the 
primary and secondary break-up processes, thus its diameter decreased with time. The variation of 
the thermal boundary layer within the droplet also changed with time. The purpose here is to 
investigate the concept of boundary layer film thickness within liquid droplets involving two 
temperatures. In the course of study, it was found that in the secondary break-up regime (i.e. the T- 
TAB regime) the droplets were so small they were heated up rapidly. Thus the current model was 
only implemented within the T-blob (i.e., the primary break-up) model, Figure 1 shows the history 
of the normalized effective film thickness within the liquid droplets. The results indicate that the 
normalized thermal boundary layer thickness is rather thick in the initial stage, and then decreases 
quickly to a small value exhibiting a reasonable physical trend. It should be noticed that the current 
model gives thinner thermal boundary layer thickness when compared with the limiting thermal 
layer thickness based on internal vortex convection model. Utilizing the model of Abramzon and 
Sirignano [8], Ra and Reitz [7] suggested that the value of thermal layer thickness be limited to 
1/2.257 of the droplet radius. It was also observed that, based on the current model, the turbulent 
diffusivity within the droplet was at least one order of magnitude greater than the laminar thermal 
diffusivity. The turbulent thermal diffusivity decreased in tandem with the decrease in the kinetic 
energy experienced by the drop. The variation of the drop surface temperature (T s ) and the bulk 
temperature of the drop (Td) for the finite conductivity (F-C) model in comparison with the bulk 
temperature calculated using the infinite conductivity (I-C) model is shown in Fig. 2. The 
temperature profiles for Td and T s were in close comparison to the results obtained by reference [7] 
for iso-octane evaporation. The variation of the normalized parent drop radius with time for a one- 



way coupled test case, with iso-octane fuel evaporating in a quiescent environment is shown in 
figure 3. The present F-C model predicted a higher rate of change in the droplet radius (i.e. shorter 
droplet life time) than the I-C model due to the higher droplet surface temperature (see Fig. 2), 
which gave a higher evaporation rate. The results of this comparison indicated that the droplet 
lifetimes predicted by the finite conductivity model are shorter than the prediction by the infinite 
conductivity model. The reason for this is due to droplet radius being evaluated based on the droplet 
surface temperature, which is higher than the interior temperature in the case of the evaporating 
droplet. The difference in rates increases with an increase in ambient temperature. The current 
model compares favourably with the models of Refs. [5] and [7], The next step in the numerical 
development of the model was to validate its efficiency for applications in practical situations. To 
this effect, the model was then incorporated into the finite volume CFD code ACE+ [14] for two- 
way coupled validation study using the Eulerian-Lagrangian methodology [15]. In the two-way 
simulation, mutual interaction between the gas phase and the evaporating spray were accounted for 
by extra source terms in the equations of change for mass and momentum. The detailed numerical 
implementation, as well as grid size and time steps issues, can be found in Balasubramanyam [13]. 

Two-Way Coupled Evaporating Spray Validation 

To evaluate the current evaporation model, the 2-D axis-symmetric subsonic low- 
evaporating spray of Yakota et al. [16] was tested for a two-way coupling CFD calculation since a 
similar non-evaporating test case was used for the T-blob/T-TAB validation study [10]. Liquid fuel 
(Tridecane, C 13 H 28 ) is injected through a single-hole nozzle into a pressurized, high temperature 
ambient N 2 environment. The test conditions for the evaporating spray are summarized in Table 1 . 
The nozzle diameter was 0.16 mm. Eqs. (1) and (2), with the same values for nozzle parameters 
except for nozzle diameter and jet velocity, were used to estimate the initial liquid turbulence 
quantities. A computational domain of 20 mm in radius and 100 mm in length discretized by a 50 
radial x 75 axial grid was used. The mesh spacing was non-uniform with refinement on the 
centreline and close to the injector. A constant time step of 2.5E-6 sec. was used with an injection 



period of 4 fi sec. The properties of liquid fuel Tridecane were taken from the NIST/JANAF 
database. Estimating the penetration of a fuel jet into an air stream is an important global property 
for model validation and is presented first. In Fig. 4, the predicted tip penetration results using the 
current finite -conductivity (F-C) evaporation model coupled with the T-blob/T-TAB atomization 
model are compared with the measured data. For reference, predictions using classical 
Blob/TAB/infinite-conductivity (I-C) model, as well as using T -blob/T -TAB/infinite-conducti vity 
model are also shown in the same figure. As indicated in [16], the infinite-conductivity evaporation 
models tend to over-predict the evaporation rate for multi-dimensional two-way coupled 
calculations, thus gave shorter tip penetration. On the other hand, the finite-conductivity model 
slows down the evaporation process, and produces larger droplets and longer penetration. It should 
be noted that the coalescence model [10] was used for all two-way coupling calculation cases. The 
coalescence model is responsible for the calculated overshoot phenomena observed in the initial 
period of injection for all simulated cases. In the downstream region, the current model gives much 
better comparison with the experimental data. The heat transfer aspects of the evaporating jet are 
shown in Fig. 5 at time of 4 fisec after injection. As can be observed, the models incorporating the 
turbulence effects in the primary and secondary atomization processes give more reasonable 
qualitative pictures when compared with the classical atomization model without liquid turbulence 
effect. Comparing Figs. 5.(b) and 5.(c), the surface temperature contours predicted by the current 
finite-conductivity model show the effect of slowing down the rate of evaporation, and qualitatively 
compared well with the more sophisticated droplet heat-up model of Bertoli and Migliaccio [17]. 

Concluding Remarks 

Based on the two-temperature film theory, a new finite-conductivity model accounting for 
droplet internal turbulence effect is developed for evaporating spray numerical calculations. The 
model is an extension of the existing T-blob/T-TAB atomization/spray model which provides the 
turbulence characteristics for estimating an effective thermal diffusivity within the droplet. Based 
on the one-way simple spray results, the model exhibits reasonable physical trends in terms of 



droplets evaporation features. The current model can be efficiently incorporated into practical spray 
combustion CFD codes. In two-way Eulerian-Lagrangian multi-dimensional full CFD simulations 
utilizing CFD-ACE+ code, the current finite-conductivity model coupled with the T-blob/T-TAB 
model shows superior performance to the conventional infinite-conductivity evaporation model, by 
comparison to evaporating spray experimental data. 
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Nomenclature 

Ad Droplet surface area (m 2 ) 

B m Gas phase mass transfer number 

Cd Discharge coefficient of injector nozzle 
(0.07 in this study) 

C p Specific heat capacity of 

mixture(KJ/kg K) 

Cfx Turbulence constant 

D Binary Diffusivity (m 2 /s) 

Kc Loss coefficient due to nozzle inlet geometry (0.45) 

K* Proportionality constant (0.23) 

k\ Turbulent kinetic energy of the liquid (m 2 /s 2 ) 

Lv Latent heat of the fuel at surface temperature (KJ/kg) 
md Droplet mass (kg) 

Nu g Gas phase Nusselt number 

Pr t Turbulent Prandtl number (0.9) 

Patm Atmospheric pressure (Pa) 

Pga S Gas phase pressure (Pa) 

qi Heat transfer rate from droplet surface to droplet core 
q g Heat transfer rate from gas phase to droplet surface 

rd Droplet radius (m), D p is diameter 

Sh g Gas phase Sherwood number 

Tb Fuel boiling temperature (K) 

Td Bulk temperature of droplet (K) 

T g Ambient gas temperature (K) 

T s Droplet surface temperature (K) 

s Area ratio at nozzle contraction, (0.0 in this study) 



U Droplet injection velocity (m/s) 

Ofeff Effective thermal diffusivity (m 2 /s) 

05am Laminar thermal diffusivity (m 2 /s) 

Gturb Turbulent thermal diffusivity (m 2 /s) 

X] Liquid thermal conductivity 

Xg Gas thermal conductivity (KW/m K) 

61 Turbulent dissipation rate (m 2 /s 3 ) 

5 e Equivalent thickness of thermal boundary layer 
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Table 1. Test Conditions for the Measurement of Yokota et al. [16] 


Case 

Pinj 

Pgas 

Tamb 

Min, 

Gas Environment 

Evap. spry 

30 Mpa 

3 Mpa 

900 K 

0.00326 kg/s 

n 2 



Figure 1 . History of effective film thickness within the evaporating 
droplet using various integration time steps. 





Temperature, 




















: : ' : V. : . 1 


mis 










Figure 5. Temperature Contours (X, Y direction in m) at 4 jisec; (a) I-C 
with Blob/TAB, (b) I-C with T-blob/T-TAB, (c) F-C with T-blob /T-TAB 





